stanEx4
options(scipen=100,digits=3)
library(cmdstanr)
options(mc.cores=parallel::detectCores())
options(cmdstanr_max_rows=1000)
library(gridExtra)execute mcmc sampling
goMCMC=function(mdl,data,smp=500,wrm=100,th=1){
mcmc=mdl$sample(
data=data,
seed=1,
chains=4,
iter_sampling=smp,
iter_warmup=wrm,
thin=th,
refresh=1000
)
mcmc
}see mcmc result and parameters
seeMCMC=function(mcmc,exc='',ch=T){ # exclude 'exc' parameters from seeing
print(mcmc)
prs=mcmc$metadata()$model_params[-1] # reject lp__
for(pr in prs){
if(grepl('^y',pr)) next # not show predictive value "y*" information
if(exc!='' && grepl(paste0('^',exc),pr)) next
drw=mcmc$draws(pr)
if(ch){
par(mfrow=c(1,4))
drw[,1,] |> plot(type='l',main='chain1',ylab=pr)
drw[,2,] |> plot(type='l',main='chain2',ylab=pr)
drw[,3,] |> plot(type='l',main='chain3',ylab=pr)
drw[,4,] |> plot(type='l',main='chain4',ylab=pr)
}
par(mfrow=c(1,2))
drw |> hist(main=pr,xlab='')
drw |> density() |> plot(main=pr)
}
}fn=function(mdl,data,smp=500,wrm=100){
mle=mdl$optimize(data=data)
print(mle)
mcmc=goMCMC(mdl,data,smp,wrm)
mcmc$metadata()$stan_variables
seeMCMC(mcmc,'m')
y0=mcmc$draws('y0')
smy0=summary(y0)
grid.arrange(
qplot(y,smy0$median,xlab='obs.',ylab='prd.')+geom_abline(intercept=0,slope=1))
par(mfrow=c(1,2))
hist(y-smy0$median,xlab='obs.-prd.',main='residual')
density(y-smy0$median) |> plot(xlab='obs.-prd.',main='residual')
m1=mcmc$draws('m1')
smm1=summary(m1)
y1=mcmc$draws('y1')
smy1=summary(y1)
xy=tibble(x=x1,m=smm1$median,ml5=smm1$q5,mu5=smm1$q95,yl5=smy1$q5,yu5=smy1$q95)
qplot(x,y,col=I('red'))+
geom_line(aes(x=x,y=ml5),xy,col='darkgray')+
geom_line(aes(x=x,y=mu5),xy,col='darkgray')+
geom_line(aes(x=x,y=yl5),xy,col='lightgray')+
geom_line(aes(x=x,y=yu5),xy,col='lightgray')+
geom_line(aes(x=x,y=m),xy,col='black')
}explanatory variable obs.x have noise
xN(x0.sx),yN(b0+b1*x0,s)
ex8-0.stan
normal regression without explanatory variable’s noise
data{
int N;
vector[N] x;
vector[N] y;
int N1;
vector[N1] x1;
}
parameters{
real b0;
real b1;
real<lower=0> s;
}
model{
y~normal(b0+b1*x,s);
}
generated quantities{
vector[N] m0;
vector[N] y0;
for(i in 1:N){
m0[i]=b0+b1*x[i];
y0[i]=normal_rng(m0[i],s);
}
vector[N1] m1;
vector[N1] y1;
for(i in 1:N1){
m1[i]=b0+b1*x1[i];
y1[i]=normal_rng(m1[i],s);
}
}ex9.stan
normal regression with explanatory variable’s noise
data{
int N;
vector[N] x;
vector[N] y;
int N1;
vector[N1] x10;
}
parameters{
real b0;
real b1;
real<lower=0> s;
real<lower=0> sx;
vector[N] x0;
}
model{
for(i in 1:N){
x[i]~normal(x0[i],sx);
y[i]~normal(b0+b1*x0[i],s);
}
}
generated quantities{
vector[N] m0;
vector[N] y0;
for(i in 1:N){
m0[i]=b0+b1*x0[i];
y0[i]=normal_rng(m0[i],s);
}
vector[N1] m1;
vector[N1] x1;
vector[N1] y1;
for(i in 1:N1){
x1[i]=normal_rng(x10[i],sx);
m1[i]=b0+b1*x10[i];
y1[i]=normal_rng(m1[i],s);
}
}n=20
x0=runif(n,0,20)
x=rnorm(n,x0,2)
y=rnorm(n,x0*2+5,2)
qplot(x,y)
n1=10
#explanatory variable do not has noise
x1=seq(min(x),max(x),length.out=n1) # new data fpr predcition
data=list(N=n,x=x,y=y,N1=n1,x1=x1)
mdl=cmdstan_model('./ex8-0.stan')
fn(mdl,data)## Initial log joint probability = -3212.4
## Iter log prob ||dx|| ||grad|| alpha alpha0 # evals Notes
## 18 -44.2348 0.000571428 0.000101815 1 1 30
## Optimization terminated normally:
## Convergence detected: relative gradient magnitude is below tolerance
## Finished in 0.2 seconds.
## variable estimate
## lp__ -44.23
## b0 6.31
## b1 1.87
## s 5.54
## m0[1] 18.23
## m0[2] 21.47
## m0[3] 13.17
## m0[4] 22.25
## m0[5] 6.33
## m0[6] 10.94
## m0[7] 43.59
## m0[8] 16.32
## m0[9] 5.80
## m0[10] 33.11
## m0[11] 26.81
## m0[12] 41.24
## m0[13] 17.51
## m0[14] 28.43
## m0[15] 28.94
## m0[16] 40.29
## m0[17] 24.08
## m0[18] 21.39
## m0[19] 15.65
## m0[20] 40.46
## y0[1] 22.51
## y0[2] 29.22
## y0[3] 14.21
## y0[4] 22.20
## y0[5] 11.34
## y0[6] 6.84
## y0[7] 46.17
## y0[8] 17.39
## y0[9] 10.52
## y0[10] 36.19
## y0[11] 28.28
## y0[12] 50.55
## y0[13] 25.44
## y0[14] 20.28
## y0[15] 30.99
## y0[16] 29.41
## y0[17] 21.73
## y0[18] 13.53
## y0[19] 6.07
## y0[20] 45.62
## m1[1] 5.80
## m1[2] 10.00
## m1[3] 14.20
## m1[4] 18.39
## m1[5] 22.59
## m1[6] 26.79
## m1[7] 30.99
## m1[8] 35.19
## m1[9] 39.39
## m1[10] 43.59
## y1[1] 8.46
## y1[2] 13.16
## y1[3] 15.73
## y1[4] 15.14
## y1[5] 19.36
## y1[6] 30.33
## y1[7] 26.02
## y1[8] 40.43
## y1[9] 43.25
## y1[10] 47.33
## Running MCMC with 4 parallel chains...
##
## Chain 1 WARNING: There aren't enough warmup iterations to fit the
## Chain 1 three stages of adaptation as currently configured.
## Chain 1 Reducing each adaptation stage to 15%/75%/10% of
## Chain 1 the given number of warmup iterations:
## Chain 1 init_buffer = 15
## Chain 1 adapt_window = 75
## Chain 1 term_buffer = 10
## Chain 1 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 1 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 1 Iteration: 600 / 600 [100%] (Sampling)
## Chain 2 WARNING: There aren't enough warmup iterations to fit the
## Chain 2 three stages of adaptation as currently configured.
## Chain 2 Reducing each adaptation stage to 15%/75%/10% of
## Chain 2 the given number of warmup iterations:
## Chain 2 init_buffer = 15
## Chain 2 adapt_window = 75
## Chain 2 term_buffer = 10
## Chain 2 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 2 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 2 Iteration: 600 / 600 [100%] (Sampling)
## Chain 3 WARNING: There aren't enough warmup iterations to fit the
## Chain 3 three stages of adaptation as currently configured.
## Chain 3 Reducing each adaptation stage to 15%/75%/10% of
## Chain 3 the given number of warmup iterations:
## Chain 3 init_buffer = 15
## Chain 3 adapt_window = 75
## Chain 3 term_buffer = 10
## Chain 3 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 3 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 3 Iteration: 600 / 600 [100%] (Sampling)
## Chain 4 WARNING: There aren't enough warmup iterations to fit the
## Chain 4 three stages of adaptation as currently configured.
## Chain 4 Reducing each adaptation stage to 15%/75%/10% of
## Chain 4 the given number of warmup iterations:
## Chain 4 init_buffer = 15
## Chain 4 adapt_window = 75
## Chain 4 term_buffer = 10
## Chain 4 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 4 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 4 Iteration: 600 / 600 [100%] (Sampling)
## Chain 1 finished in 0.1 seconds.
## Chain 2 finished in 0.1 seconds.
## Chain 3 finished in 0.0 seconds.
## Chain 4 finished in 0.0 seconds.
##
## All 4 chains finished successfully.
## Mean chain execution time: 0.1 seconds.
## Total execution time: 0.4 seconds.
##
## variable mean median sd mad q5 q95 rhat ess_bulk ess_tail
## lp__ -44.18 -43.76 1.44 1.08 -47.10 -42.70 1.00 405 425
## b0 6.70 6.56 2.71 2.19 2.40 11.36 1.01 265 204
## b1 1.83 1.83 0.24 0.22 1.43 2.24 1.01 357 409
## s 6.32 6.16 1.20 1.05 4.69 8.53 1.00 951 802
## m0[1] 18.39 18.40 1.65 1.45 15.68 21.24 1.01 336 335
## m0[2] 21.58 21.62 1.50 1.41 19.15 24.00 1.01 458 518
## m0[3] 13.43 13.38 2.03 1.73 10.16 16.97 1.01 276 225
## m0[4] 22.34 22.39 1.49 1.40 20.00 24.73 1.01 511 590
## m0[5] 6.72 6.58 2.71 2.19 2.42 11.39 1.01 265 204
## m0[6] 11.24 11.15 2.24 1.87 7.62 15.08 1.01 269 217
## m0[7] 43.27 43.27 2.88 2.80 38.51 47.97 1.00 999 933
## m0[8] 16.52 16.50 1.78 1.52 13.72 19.58 1.01 303 256
## m0[9] 6.20 6.06 2.77 2.24 1.82 10.95 1.01 264 204
## m0[10] 33.00 32.98 1.86 1.80 29.88 36.00 1.00 2059 1320
## m0[11] 26.81 26.83 1.50 1.45 24.30 29.24 1.00 1240 1072
## m0[12] 40.97 40.95 2.63 2.55 36.59 45.27 1.00 1190 1114
## m0[13] 17.69 17.68 1.69 1.47 14.93 20.66 1.01 321 298
## m0[14] 28.40 28.41 1.56 1.51 25.77 30.90 1.00 1780 1190
## m0[15] 28.91 28.92 1.59 1.52 26.24 31.44 1.00 1885 1207
## m0[16] 40.03 40.02 2.53 2.47 35.81 44.16 1.00 1289 1058
## m0[17] 24.13 24.15 1.46 1.41 21.75 26.54 1.00 698 754
## m0[18] 21.50 21.54 1.51 1.40 19.06 23.93 1.01 452 518
## m0[19] 15.87 15.84 1.83 1.56 12.99 19.01 1.01 295 244
## m0[20] 40.20 40.18 2.55 2.48 35.96 44.37 1.00 1269 1058
## y0[1] 18.28 18.31 6.72 6.38 7.25 29.63 1.00 1367 1751
## y0[2] 21.51 21.33 6.60 6.36 10.71 32.22 1.00 2103 1946
## y0[3] 13.47 13.66 6.72 6.38 2.69 24.19 1.00 1377 1791
## y0[4] 22.24 22.16 6.65 6.52 11.69 33.54 1.00 1465 1581
## y0[5] 6.69 6.83 6.96 6.82 -4.66 18.16 1.01 976 1169
## y0[6] 11.17 11.31 6.82 6.51 0.23 22.15 1.00 1277 934
## y0[7] 43.55 43.44 7.18 6.42 32.38 55.60 1.00 1714 1791
## y0[8] 16.29 16.36 6.66 6.23 5.28 27.30 1.00 1279 1661
## y0[9] 6.21 6.17 6.94 6.47 -5.33 17.54 1.00 875 994
## y0[10] 32.90 32.88 6.76 6.53 21.66 44.03 1.00 1898 1930
## y0[11] 26.87 26.91 6.71 6.57 16.33 37.58 1.00 2004 1974
## y0[12] 41.23 41.02 6.99 6.63 30.25 52.79 1.00 1955 1856
## y0[13] 17.60 17.79 6.86 6.67 6.45 28.43 1.01 1222 1330
## y0[14] 28.49 28.47 6.68 6.32 17.66 39.27 1.00 2124 1970
## y0[15] 28.80 28.59 6.75 6.23 17.62 40.05 1.00 2166 1907
## y0[16] 40.05 39.94 7.04 6.85 28.60 51.76 1.00 2077 1955
## y0[17] 24.04 24.01 6.58 6.21 13.28 34.51 1.00 1736 1863
## y0[18] 21.63 21.61 6.68 6.49 10.88 32.62 1.00 1783 1876
## y0[19] 15.79 16.00 6.59 6.08 4.72 26.17 1.00 1529 1731
## y0[20] 40.22 40.23 6.88 6.46 28.71 51.18 1.00 1985 1838
## m1[1] 6.20 6.06 2.77 2.24 1.82 10.95 1.01 264 204
## m1[2] 10.32 10.23 2.33 1.93 6.60 14.33 1.01 267 215
## m1[3] 14.44 14.41 1.94 1.65 11.31 17.85 1.01 282 236
## m1[4] 18.56 18.56 1.64 1.45 15.86 21.37 1.01 340 335
## m1[5] 22.68 22.72 1.48 1.40 20.35 25.08 1.01 539 612
## m1[6] 26.79 26.82 1.50 1.45 24.29 29.23 1.00 1235 1072
## m1[7] 30.91 30.91 1.71 1.63 28.06 33.71 1.00 2088 1297
## m1[8] 35.03 35.02 2.03 1.96 31.62 38.36 1.00 1845 1286
## m1[9] 39.15 39.14 2.44 2.36 35.07 43.10 1.00 1402 1037
## m1[10] 43.27 43.27 2.88 2.80 38.51 47.97 1.00 999 933
## y1[1] 6.25 6.31 6.84 6.60 -5.03 17.77 1.00 1070 1388
## y1[2] 10.18 10.37 6.81 6.66 -1.02 21.67 1.00 1019 1279
## y1[3] 14.50 14.47 6.77 6.52 3.63 25.58 1.00 1171 1522
## y1[4] 18.65 18.57 6.82 6.49 7.73 29.42 1.00 1836 1923
## y1[5] 22.57 22.71 6.85 6.44 11.53 33.82 1.00 1863 1933
## y1[6] 26.86 26.92 6.49 6.39 16.25 37.07 1.00 2135 2009
## y1[7] 30.89 30.80 6.72 6.35 20.04 41.83 1.00 1600 1907
## y1[8] 34.92 34.99 6.91 6.41 23.67 46.13 1.00 1970 1894
## y1[9] 39.28 39.23 6.77 6.34 27.99 50.30 1.00 1933 1953
## y1[10] 43.37 43.53 7.14 6.89 31.65 54.73 1.00 1618 1448
#explanatory variable has noise
x10=seq(min(x),max(x),length.out=n1) # new data fpr predcition
data=list(N=n,x=x,y=y,N1=n1,x10=x10)
mdl=cmdstan_model('./ex9.stan')
mcmc=goMCMC(mdl,data,wrm=500,smp=1000)## Running MCMC with 4 parallel chains...
##
## Chain 1 Iteration: 1 / 1500 [ 0%] (Warmup)
## Chain 2 Iteration: 1 / 1500 [ 0%] (Warmup)
## Chain 3 Iteration: 1 / 1500 [ 0%] (Warmup)
## Chain 4 Iteration: 1 / 1500 [ 0%] (Warmup)
## Chain 1 Iteration: 501 / 1500 [ 33%] (Sampling)
## Chain 3 Iteration: 501 / 1500 [ 33%] (Sampling)
## Chain 4 Iteration: 501 / 1500 [ 33%] (Sampling)
## Chain 2 Iteration: 501 / 1500 [ 33%] (Sampling)
## Chain 1 Iteration: 1500 / 1500 [100%] (Sampling)
## Chain 3 Iteration: 1500 / 1500 [100%] (Sampling)
## Chain 1 finished in 0.5 seconds.
## Chain 3 finished in 0.5 seconds.
## Chain 2 Iteration: 1500 / 1500 [100%] (Sampling)
## Chain 2 finished in 0.6 seconds.
## Chain 4 Iteration: 1500 / 1500 [100%] (Sampling)
## Chain 4 finished in 0.6 seconds.
##
## All 4 chains finished successfully.
## Mean chain execution time: 0.6 seconds.
## Total execution time: 0.7 seconds.seeMCMC(mcmc,'[m,x]')## variable mean median sd mad q5 q95 rhat ess_bulk ess_tail
## lp__ -57.13 -58.46 10.96 9.52 -73.12 -36.69 1.02 147 159
## b0 6.73 6.67 2.86 2.77 2.26 11.46 1.01 873 1170
## b1 1.82 1.82 0.26 0.24 1.40 2.24 1.01 728 979
## s 4.36 4.42 1.93 2.11 1.24 7.49 1.06 65 171
## sx 2.36 2.37 1.27 1.29 0.42 4.47 1.04 109 148
## x0[1] 4.98 5.22 1.87 1.83 1.87 7.59 1.02 326 1386
## x0[2] 9.42 9.32 1.72 1.76 6.82 12.30 1.02 384 1902
## x0[3] 1.94 2.28 2.11 2.06 -1.66 4.75 1.02 286 883
## x0[4] 8.74 8.69 1.54 1.19 6.21 11.34 1.01 3254 1758
## x0[5] 1.00 0.93 1.80 1.67 -1.76 3.89 1.01 846 1874
## x0[6] 2.63 2.63 1.67 1.38 -0.14 5.31 1.01 1881 1429
## x0[7] 18.38 18.54 1.96 2.08 15.11 21.20 1.01 511 1670
## x0[8] 5.62 5.59 1.54 1.26 3.22 8.04 1.00 2340 1800
## x0[9] -0.59 -0.42 1.75 1.41 -3.71 2.05 1.00 1395 907
## x0[10] 14.19 14.20 1.56 1.33 11.68 16.82 1.01 2232 1896
## x0[11] 14.01 13.83 2.63 3.18 10.51 18.30 1.04 146 1185
## x0[12] 18.40 18.43 1.66 1.42 15.77 21.07 1.00 1382 1201
## x0[13] 1.81 2.13 3.39 4.15 -3.84 6.26 1.04 121 692
## x0[14] 13.38 13.17 1.89 1.91 10.70 16.65 1.02 247 1281
## x0[15] 12.49 12.40 1.54 1.25 10.10 15.12 1.00 1669 1386
## x0[16] 18.49 18.30 1.71 1.40 16.00 21.46 1.01 1338 1488
## x0[17] 10.32 10.20 1.64 1.42 7.90 13.01 1.01 900 2024
## x0[18] 9.92 9.83 1.97 2.20 7.05 13.16 1.03 146 1397
## x0[19] 3.62 3.83 1.88 1.79 0.38 6.25 1.02 305 1202
## x0[20] 18.74 18.54 1.85 1.41 16.07 21.85 1.01 1256 1069
## m0[1] 15.93 15.80 3.16 3.29 11.08 21.04 1.01 507 2554
## m0[2] 23.85 23.88 3.09 3.27 18.84 28.70 1.02 338 2115
## m0[3] 10.43 10.27 3.56 3.70 5.02 16.32 1.02 377 2073
## m0[4] 22.64 22.66 2.72 2.33 18.14 27.02 1.01 4145 2190
## m0[5] 8.68 8.88 3.38 3.41 3.03 13.77 1.01 585 2142
## m0[6] 11.64 11.57 2.96 2.62 6.84 16.56 1.01 2756 1891
## m0[7] 40.06 39.83 3.74 3.92 34.63 46.38 1.02 291 1585
## m0[8] 17.02 17.05 2.75 2.42 12.57 21.48 1.00 2800 2598
## m0[9] 5.84 5.71 3.11 2.74 0.86 11.10 1.00 2283 1892
## m0[10] 32.49 32.45 2.81 2.46 27.93 37.22 1.00 3554 2242
## m0[11] 32.09 31.99 4.50 5.66 25.17 38.85 1.03 113 1931
## m0[12] 40.07 40.01 2.99 2.65 35.18 45.05 1.00 1719 1999
## m0[13] 10.23 10.22 5.74 7.51 1.82 18.94 1.04 128 1210
## m0[14] 30.98 31.05 3.22 3.50 25.76 36.05 1.02 308 2435
## m0[15] 29.41 29.46 2.72 2.43 24.92 33.67 1.00 3499 2413
## m0[16] 40.22 40.25 2.94 2.73 35.33 45.08 1.00 3383 2444
## m0[17] 25.48 25.56 2.88 2.70 20.81 30.07 1.01 1079 2531
## m0[18] 24.75 24.75 3.52 3.88 19.20 30.23 1.03 159 2588
## m0[19] 13.45 13.28 3.23 3.22 8.46 18.74 1.01 546 2091
## m0[20] 40.66 40.75 3.13 2.79 35.26 45.71 1.00 3493 2437
## y0[1] 15.85 15.14 5.75 4.86 7.33 26.20 1.01 1876 2915
## y0[2] 24.02 24.53 5.71 4.99 13.90 32.78 1.01 1560 2453
## y0[3] 10.44 9.68 5.81 5.21 2.02 20.76 1.01 1236 3227
## y0[4] 22.49 22.63 5.44 4.61 13.40 31.26 1.01 3869 2910
## y0[5] 8.64 9.08 5.99 5.07 -2.08 17.55 1.01 1567 2687
## y0[6] 11.62 11.63 5.60 4.74 2.19 20.64 1.01 4354 3407
## y0[7] 39.93 39.16 6.10 5.39 31.26 50.63 1.01 1177 1901
## y0[8] 17.09 17.15 5.39 4.49 8.10 26.07 1.01 3307 3332
## y0[9] 5.84 5.74 5.66 4.60 -3.47 15.43 1.01 3313 2505
## y0[10] 32.59 32.44 5.53 4.48 23.79 42.04 1.01 3861 3171
## y0[11] 31.92 32.80 6.55 6.44 20.11 40.84 1.02 289 2071
## y0[12] 40.09 39.79 5.60 4.51 30.91 49.79 1.01 2897 3274
## y0[13] 10.25 9.20 7.55 7.87 0.36 23.75 1.02 264 1353
## y0[14] 30.84 31.65 5.78 4.94 20.55 39.18 1.01 1273 2721
## y0[15] 29.39 29.58 5.48 4.44 20.04 38.04 1.01 4142 2891
## y0[16] 40.18 40.20 5.51 4.63 31.23 49.36 1.01 3823 3145
## y0[17] 25.60 26.03 5.52 4.66 15.94 33.97 1.01 2616 2582
## y0[18] 24.65 25.49 5.96 5.31 13.92 33.43 1.01 986 2472
## y0[19] 13.35 12.70 5.74 4.82 4.92 23.58 1.01 1807 2595
## y0[20] 40.62 40.91 5.80 4.79 30.85 49.71 1.01 3553 2938
## m1[1] 6.23 6.17 2.92 2.82 1.67 11.07 1.01 866 1135
## m1[2] 10.33 10.29 2.44 2.33 6.44 14.34 1.01 936 1265
## m1[3] 14.42 14.38 2.01 1.92 11.15 17.79 1.00 1080 1762
## m1[4] 18.52 18.51 1.66 1.58 15.79 21.23 1.00 1356 2225
## m1[5] 22.61 22.61 1.47 1.39 20.18 24.96 1.00 1715 2550
## m1[6] 26.70 26.70 1.50 1.44 24.26 29.08 1.00 1658 2398
## m1[7] 30.80 30.78 1.72 1.66 27.93 33.56 1.00 1229 1923
## m1[8] 34.89 34.87 2.09 2.04 31.47 38.24 1.01 975 1417
## m1[9] 38.99 38.96 2.54 2.44 34.84 43.11 1.01 866 1304
## m1[10] 43.08 43.03 3.02 2.87 38.15 48.07 1.01 814 1295
## x1[1] -0.27 -0.26 2.71 1.89 -4.68 4.12 1.01 3822 3082
## x1[2] 2.01 1.98 2.65 1.82 -2.10 6.32 1.02 3924 2630
## x1[3] 4.23 4.24 2.66 1.91 -0.24 8.72 1.02 4081 1948
## x1[4] 6.42 6.46 2.65 1.88 2.01 10.71 1.01 3486 2368
## x1[5] 8.71 8.70 2.59 1.83 4.46 13.06 1.01 4154 2681
## x1[6] 10.94 10.97 2.72 1.95 6.42 15.17 1.02 4062 2458
## x1[7] 13.23 13.22 2.63 1.96 8.91 17.63 1.01 4086 2501
## x1[8] 15.45 15.46 2.66 1.94 11.05 19.77 1.01 3940 2459
## x1[9] 17.71 17.72 2.71 1.95 13.39 22.06 1.01 4017 2355
## x1[10] 19.94 19.94 2.65 1.86 15.59 24.18 1.01 4498 2475
## y1[1] 6.24 6.24 5.67 4.79 -2.96 15.20 1.01 1932 2447
## y1[2] 10.46 10.49 5.25 4.46 1.83 19.13 1.01 2863 2619
## y1[3] 14.48 14.51 5.14 4.31 6.07 23.02 1.01 3040 2695
## y1[4] 18.38 18.52 5.13 4.20 9.80 26.61 1.01 2842 2651
## y1[5] 22.73 22.74 5.00 4.04 14.65 31.07 1.01 3721 3043
## y1[6] 26.60 26.59 4.99 3.98 18.45 34.79 1.01 3139 2788
## y1[7] 30.70 30.70 4.96 4.22 22.34 38.87 1.01 3151 2833
## y1[8] 34.96 34.86 5.32 4.48 26.42 44.01 1.01 2489 2777
## y1[9] 38.99 38.97 5.26 4.47 30.32 47.55 1.01 2444 3103
## y1[10] 43.16 43.23 5.72 5.00 33.58 52.59 1.01 1520 2616
y0=mcmc$draws('y0')
smy0=summary(y0)
grid.arrange(
qplot(y,smy0$median,xlab='obs.',ylab='prd.')+geom_abline(intercept=0,slope=1))
par(mfrow=c(1,2))
hist(y-smy0$median,xlab='obs.-prd.',main='residual')
density(y-smy0$median) |> plot(xlab='obs.-prd.',main='residual')
m1=mcmc$draws('m1')
smm1=summary(m1)
x1=mcmc$draws('x1')
smx1=summary(x1)
y1=mcmc$draws('y1')
smy1=summary(y1)
xy=tibble(x=smx1$median,m=smm1$median,ml5=smm1$q5,mu5=smm1$q95,yl5=smy1$q5,yu5=smy1$q95)
qplot(x,y,col=I('red'))+
geom_line(aes(x=x,y=ml5),xy,col='darkgray')+
geom_line(aes(x=x,y=mu5),xy,col='darkgray')+
geom_line(aes(x=x,y=yl5),xy,col='lightgray')+
geom_line(aes(x=x,y=yu5),xy,col='lightgray')+
geom_line(aes(x=x,y=m),xy,col='black')
outlier
ex10.stan
objective variable have outlier, y~cauchy(b0+b1*x,s)
data{
int N;
vector[N] x;
vector[N] y;
int N1;
vector[N1] x1;
}
parameters{
real b0;
real b1;
real<lower=0> s;
}
model{
y~cauchy(b0+b1*x,s);
}
generated quantities{
vector[N] m0;
vector[N] y0;
for(i in 1:N){
m0[i]=b0+b1*x[i];
y0[i]=cauchy_rng(m0[i],s);
}
vector[N1] m1;
vector[N1] y1;
for(i in 1:N1){
m1[i]=b0+b1*x1[i];
y1[i]=cauchy_rng(m1[i],s);
}
}n=20
x=runif(n,0,9)
y=rnorm(n,x*2+5,1)
x[1]=3
y[1]=25
qplot(x,y)
n1=10
x1=seq(min(x),max(x),length.out=n1) # new data fpr predcition
data=list(N=n,x=x,y=y,N1=n1,x1=x1)
mdl=cmdstan_model('./ex8-0.stan')
fn(mdl,data)## Initial log joint probability = -2635.41
## Iter log prob ||dx|| ||grad|| alpha alpha0 # evals Notes
## 16 -33.4085 0.000232464 0.000110767 0.3378 0.9899 26
## Optimization terminated normally:
## Convergence detected: relative gradient magnitude is below tolerance
## Finished in 0.1 seconds.
## variable estimate
## lp__ -33.41
## b0 6.76
## b1 1.80
## s 3.22
## m0[1] 12.15
## m0[2] 22.66
## m0[3] 16.16
## m0[4] 17.89
## m0[5] 16.68
## m0[6] 21.87
## m0[7] 22.46
## m0[8] 18.49
## m0[9] 15.72
## m0[10] 8.08
## m0[11] 15.38
## m0[12] 17.28
## m0[13] 20.49
## m0[14] 9.05
## m0[15] 12.83
## m0[16] 15.95
## m0[17] 19.81
## m0[18] 21.84
## m0[19] 9.24
## m0[20] 21.01
## y0[1] 12.49
## y0[2] 21.55
## y0[3] 16.14
## y0[4] 13.73
## y0[5] 18.04
## y0[6] 23.37
## y0[7] 22.50
## y0[8] 19.11
## y0[9] 17.11
## y0[10] 8.63
## y0[11] 11.79
## y0[12] 18.73
## y0[13] 17.70
## y0[14] 5.52
## y0[15] 13.27
## y0[16] 15.39
## y0[17] 12.43
## y0[18] 18.51
## y0[19] 11.25
## y0[20] 17.72
## m1[1] 8.08
## m1[2] 9.70
## m1[3] 11.32
## m1[4] 12.94
## m1[5] 14.56
## m1[6] 16.18
## m1[7] 17.80
## m1[8] 19.42
## m1[9] 21.04
## m1[10] 22.66
## y1[1] 7.39
## y1[2] 9.52
## y1[3] 14.49
## y1[4] 15.14
## y1[5] 12.61
## y1[6] 18.93
## y1[7] 8.08
## y1[8] 22.46
## y1[9] 22.23
## y1[10] 20.24
## Running MCMC with 4 parallel chains...
##
## Chain 1 WARNING: There aren't enough warmup iterations to fit the
## Chain 1 three stages of adaptation as currently configured.
## Chain 1 Reducing each adaptation stage to 15%/75%/10% of
## Chain 1 the given number of warmup iterations:
## Chain 1 init_buffer = 15
## Chain 1 adapt_window = 75
## Chain 1 term_buffer = 10
## Chain 1 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 1 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 1 Iteration: 600 / 600 [100%] (Sampling)
## Chain 2 WARNING: There aren't enough warmup iterations to fit the
## Chain 2 three stages of adaptation as currently configured.
## Chain 2 Reducing each adaptation stage to 15%/75%/10% of
## Chain 2 the given number of warmup iterations:
## Chain 2 init_buffer = 15
## Chain 2 adapt_window = 75
## Chain 2 term_buffer = 10
## Chain 2 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 2 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 2 Iteration: 600 / 600 [100%] (Sampling)
## Chain 3 WARNING: There aren't enough warmup iterations to fit the
## Chain 3 three stages of adaptation as currently configured.
## Chain 3 Reducing each adaptation stage to 15%/75%/10% of
## Chain 3 the given number of warmup iterations:
## Chain 3 init_buffer = 15
## Chain 3 adapt_window = 75
## Chain 3 term_buffer = 10
## Chain 3 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 3 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 3 Iteration: 600 / 600 [100%] (Sampling)
## Chain 4 WARNING: There aren't enough warmup iterations to fit the
## Chain 4 three stages of adaptation as currently configured.
## Chain 4 Reducing each adaptation stage to 15%/75%/10% of
## Chain 4 the given number of warmup iterations:
## Chain 4 init_buffer = 15
## Chain 4 adapt_window = 75
## Chain 4 term_buffer = 10
## Chain 4 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 4 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 4 Iteration: 600 / 600 [100%] (Sampling)
## Chain 1 finished in 0.0 seconds.
## Chain 2 finished in 0.0 seconds.
## Chain 3 finished in 0.1 seconds.
## Chain 4 finished in 0.1 seconds.
##
## All 4 chains finished successfully.
## Mean chain execution time: 0.0 seconds.
## Total execution time: 0.2 seconds.
##
## variable mean median sd mad q5 q95 rhat ess_bulk ess_tail
## lp__ -33.91 -33.59 1.33 1.14 -36.51 -32.41 1.00 550 777
## b0 6.70 6.70 2.09 2.06 3.27 9.97 1.01 253 326
## b1 1.81 1.80 0.34 0.34 1.25 2.36 1.01 300 418
## s 3.69 3.58 0.67 0.60 2.77 4.94 1.00 1290 1290
## m0[1] 12.12 12.11 1.23 1.20 10.13 14.11 1.01 286 379
## m0[2] 22.69 22.69 1.40 1.35 20.43 25.04 1.00 1044 1332
## m0[3] 16.16 16.17 0.87 0.85 14.74 17.62 1.00 801 1018
## m0[4] 17.89 17.92 0.89 0.86 16.44 19.34 1.00 1869 1569
## m0[5] 16.68 16.70 0.86 0.84 15.24 18.12 1.00 1033 1350
## m0[6] 21.89 21.90 1.29 1.24 19.82 24.04 1.00 1256 1423
## m0[7] 22.49 22.50 1.37 1.32 20.28 24.81 1.00 1088 1332
## m0[8] 18.50 18.52 0.92 0.90 17.02 20.01 1.00 2244 1623
## m0[9] 15.71 15.72 0.89 0.86 14.27 17.18 1.01 652 987
## m0[10] 8.03 8.06 1.86 1.83 4.99 10.91 1.01 253 329
## m0[11] 15.37 15.38 0.90 0.87 13.92 16.91 1.01 563 817
## m0[12] 17.28 17.30 0.87 0.83 15.82 18.70 1.00 1396 1415
## m0[13] 20.51 20.52 1.11 1.07 18.69 22.32 1.00 1838 1545
## m0[14] 9.00 9.02 1.70 1.67 6.21 11.65 1.01 255 314
## m0[15] 12.81 12.78 1.14 1.13 10.96 14.68 1.01 307 511
## m0[16] 15.94 15.95 0.88 0.85 14.52 17.41 1.00 729 1029
## m0[17] 19.83 19.85 1.03 0.99 18.15 21.50 1.00 2124 1556
## m0[18] 21.87 21.87 1.29 1.24 19.80 24.01 1.00 1263 1444
## m0[19] 9.20 9.22 1.67 1.63 6.45 11.79 1.01 256 314
## m0[20] 21.03 21.03 1.17 1.13 19.11 22.97 1.00 1609 1514
## y0[1] 12.10 11.99 4.01 3.85 5.63 18.87 1.00 1336 1448
## y0[2] 22.59 22.71 4.03 3.91 15.84 28.89 1.00 1987 1845
## y0[3] 16.10 16.08 3.82 3.70 9.92 22.39 1.00 1727 1852
## y0[4] 17.86 17.90 3.97 3.87 11.46 24.29 1.00 2002 1995
## y0[5] 16.75 16.69 3.92 3.69 10.61 23.27 1.00 1773 1895
## y0[6] 21.98 21.95 4.17 3.95 15.24 28.79 1.00 1692 1592
## y0[7] 22.53 22.48 3.95 3.90 16.17 29.34 1.00 1917 1778
## y0[8] 18.53 18.53 3.86 3.79 12.17 24.77 1.00 1918 2053
## y0[9] 15.65 15.52 3.82 3.73 9.44 21.88 1.00 2156 1713
## y0[10] 8.14 8.13 4.17 4.04 1.39 15.20 1.00 1115 1330
## y0[11] 15.17 15.24 3.81 3.73 8.86 21.35 1.00 1730 1972
## y0[12] 17.46 17.44 3.85 3.72 11.30 23.98 1.00 1990 1912
## y0[13] 20.54 20.50 3.94 3.73 14.35 27.27 1.00 2021 1725
## y0[14] 8.94 8.99 3.99 4.02 2.36 15.25 1.00 835 1614
## y0[15] 12.78 12.81 3.95 3.79 6.30 18.98 1.00 1243 1744
## y0[16] 15.84 15.88 3.71 3.59 9.89 22.04 1.00 1991 1681
## y0[17] 20.01 20.00 3.88 3.70 13.65 26.25 1.00 2029 1916
## y0[18] 21.94 21.92 3.94 3.72 15.82 28.36 1.00 2009 1934
## y0[19] 9.19 9.18 4.07 4.02 2.36 15.87 1.00 894 1333
## y0[20] 20.99 20.85 4.00 3.85 14.76 27.58 1.00 1831 1940
## m1[1] 8.03 8.06 1.86 1.83 4.99 10.91 1.01 253 329
## m1[2] 9.66 9.68 1.60 1.56 7.04 12.15 1.01 258 341
## m1[3] 11.29 11.28 1.35 1.32 9.10 13.46 1.01 271 361
## m1[4] 12.92 12.89 1.13 1.12 11.10 14.77 1.01 312 501
## m1[5] 14.55 14.54 0.96 0.91 12.97 16.16 1.01 424 640
## m1[6] 16.17 16.19 0.87 0.85 14.76 17.64 1.00 807 1075
## m1[7] 17.80 17.82 0.88 0.85 16.35 19.25 1.00 1790 1568
## m1[8] 19.43 19.45 1.00 0.96 17.85 21.06 1.00 2248 1543
## m1[9] 21.06 21.06 1.18 1.13 19.14 23.00 1.00 1597 1514
## m1[10] 22.69 22.69 1.40 1.35 20.43 25.04 1.00 1044 1332
## y1[1] 7.99 8.06 4.15 3.91 1.35 14.76 1.00 1000 1345
## y1[2] 9.62 9.60 4.07 3.86 3.01 16.46 1.00 1122 1682
## y1[3] 11.34 11.31 4.02 3.99 4.81 17.89 1.00 1308 1694
## y1[4] 12.93 12.83 4.03 3.75 6.52 19.50 1.00 1607 1841
## y1[5] 14.47 14.45 3.87 3.76 7.96 20.93 1.00 1727 2012
## y1[6] 16.21 16.25 3.84 3.56 10.01 22.41 1.00 1884 1800
## y1[7] 17.68 17.58 3.82 3.79 11.38 24.07 1.00 1964 1754
## y1[8] 19.45 19.60 3.90 3.77 13.09 25.65 1.00 2213 1971
## y1[9] 21.07 20.99 3.97 3.62 14.46 27.74 1.00 1983 1775
## y1[10] 22.67 22.63 4.01 3.91 16.22 29.22 1.00 1972 1971
mdl=cmdstan_model('./ex10.stan')
fn(mdl,data)## Initial log joint probability = -72.2134
## Iter log prob ||dx|| ||grad|| alpha alpha0 # evals Notes
## 23 -15.9805 0.000356175 0.00111771 1 1 34
## Optimization terminated normally:
## Convergence detected: relative gradient magnitude is below tolerance
## Finished in 0.1 seconds.
## variable estimate
## lp__ -15.98
## b0 4.64
## b1 2.10
## s 0.68
## m0[1] 10.94
## m0[2] 23.23
## m0[3] 15.64
## m0[4] 17.65
## m0[5] 16.24
## m0[6] 22.30
## m0[7] 22.99
## m0[8] 18.36
## m0[9] 15.11
## m0[10] 6.19
## m0[11] 14.72
## m0[12] 16.94
## m0[13] 20.69
## m0[14] 7.32
## m0[15] 11.74
## m0[16] 15.38
## m0[17] 19.90
## m0[18] 22.27
## m0[19] 7.55
## m0[20] 21.30
## y0[1] 11.39
## y0[2] 22.40
## y0[3] 16.49
## y0[4] 17.86
## y0[5] -1.41
## y0[6] 23.63
## y0[7] 23.74
## y0[8] 17.73
## y0[9] -9.13
## y0[10] 3.65
## y0[11] 14.93
## y0[12] 17.02
## y0[13] 21.37
## y0[14] 7.58
## y0[15] 12.74
## y0[16] 14.20
## y0[17] 19.94
## y0[18] 22.96
## y0[19] 0.50
## y0[20] 21.10
## m1[1] 6.19
## m1[2] 8.08
## m1[3] 9.97
## m1[4] 11.87
## m1[5] 13.76
## m1[6] 15.65
## m1[7] 17.55
## m1[8] 19.44
## m1[9] 21.33
## m1[10] 23.23
## y1[1] 3.82
## y1[2] 6.63
## y1[3] 9.40
## y1[4] 11.12
## y1[5] 13.01
## y1[6] 17.35
## y1[7] 18.27
## y1[8] 16.62
## y1[9] 20.10
## y1[10] 23.67
## Running MCMC with 4 parallel chains...
##
## Chain 1 WARNING: There aren't enough warmup iterations to fit the
## Chain 1 three stages of adaptation as currently configured.
## Chain 1 Reducing each adaptation stage to 15%/75%/10% of
## Chain 1 the given number of warmup iterations:
## Chain 1 init_buffer = 15
## Chain 1 adapt_window = 75
## Chain 1 term_buffer = 10
## Chain 1 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 1 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 1 Iteration: 600 / 600 [100%] (Sampling)
## Chain 2 WARNING: There aren't enough warmup iterations to fit the
## Chain 2 three stages of adaptation as currently configured.
## Chain 2 Reducing each adaptation stage to 15%/75%/10% of
## Chain 2 the given number of warmup iterations:
## Chain 2 init_buffer = 15
## Chain 2 adapt_window = 75
## Chain 2 term_buffer = 10
## Chain 2 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 2 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 2 Iteration: 600 / 600 [100%] (Sampling)
## Chain 3 WARNING: There aren't enough warmup iterations to fit the
## Chain 3 three stages of adaptation as currently configured.
## Chain 3 Reducing each adaptation stage to 15%/75%/10% of
## Chain 3 the given number of warmup iterations:
## Chain 3 init_buffer = 15
## Chain 3 adapt_window = 75
## Chain 3 term_buffer = 10
## Chain 3 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 3 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 3 Iteration: 600 / 600 [100%] (Sampling)
## Chain 4 WARNING: There aren't enough warmup iterations to fit the
## Chain 4 three stages of adaptation as currently configured.
## Chain 4 Reducing each adaptation stage to 15%/75%/10% of
## Chain 4 the given number of warmup iterations:
## Chain 4 init_buffer = 15
## Chain 4 adapt_window = 75
## Chain 4 term_buffer = 10
## Chain 4 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 4 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 4 Iteration: 600 / 600 [100%] (Sampling)
## Chain 1 finished in 0.1 seconds.
## Chain 2 finished in 0.1 seconds.
## Chain 3 finished in 0.1 seconds.
## Chain 4 finished in 0.1 seconds.
##
## All 4 chains finished successfully.
## Mean chain execution time: 0.1 seconds.
## Total execution time: 0.2 seconds.
##
## variable mean median sd mad q5 q95 rhat ess_bulk ess_tail
## lp__ -18.00 -17.60 1.42 1.07 -20.92 -16.51 1.01 589 517
## b0 4.62 4.61 0.61 0.45 3.63 5.57 1.01 488 366
## b1 2.10 2.10 0.11 0.10 1.92 2.27 1.01 646 477
## s 0.85 0.81 0.27 0.27 0.47 1.34 1.00 607 666
## m0[1] 10.90 10.91 0.35 0.29 10.33 11.44 1.01 583 450
## m0[2] 23.15 23.16 0.55 0.57 22.24 24.01 1.00 1407 1387
## m0[3] 15.58 15.60 0.30 0.30 15.05 16.03 1.00 1577 974
## m0[4] 17.59 17.61 0.34 0.34 17.01 18.12 1.00 2101 1464
## m0[5] 16.18 16.20 0.31 0.31 15.64 16.65 1.00 1823 1068
## m0[6] 22.22 22.24 0.51 0.53 21.38 23.02 1.00 1546 1403
## m0[7] 22.91 22.93 0.54 0.56 22.02 23.76 1.00 1440 1450
## m0[8] 18.30 18.31 0.36 0.36 17.68 18.85 1.00 2092 1476
## m0[9] 15.06 15.08 0.30 0.29 14.53 15.50 1.00 1379 857
## m0[10] 6.16 6.16 0.53 0.40 5.29 7.00 1.01 487 349
## m0[11] 14.67 14.69 0.30 0.29 14.15 15.12 1.00 1244 747
## m0[12] 16.89 16.91 0.32 0.32 16.32 17.38 1.00 2026 1527
## m0[13] 20.62 20.63 0.44 0.46 19.88 21.32 1.00 1832 1564
## m0[14] 7.29 7.29 0.49 0.37 6.49 8.05 1.01 492 369
## m0[15] 11.70 11.71 0.33 0.28 11.15 12.19 1.01 639 452
## m0[16] 15.33 15.35 0.30 0.30 14.81 15.78 1.00 1479 962
## m0[17] 19.83 19.85 0.41 0.42 19.12 20.47 1.00 1969 1675
## m0[18] 22.20 22.22 0.51 0.53 21.36 22.99 1.00 1551 1468
## m0[19] 7.51 7.52 0.48 0.36 6.75 8.26 1.01 494 369
## m0[20] 21.23 21.24 0.47 0.48 20.44 21.96 1.00 1718 1617
## y0[1] 11.68 10.89 99.38 1.32 5.91 15.95 1.00 1884 1638
## y0[2] 22.68 23.17 30.39 1.40 18.04 28.52 1.00 2151 1968
## y0[3] 15.95 15.59 22.87 1.21 10.46 22.26 1.00 1873 1959
## y0[4] 17.36 17.62 19.24 1.31 12.63 22.81 1.00 2060 1854
## y0[5] 16.54 16.16 78.53 1.36 10.38 22.05 1.00 1765 1924
## y0[6] 30.65 22.25 404.41 1.38 17.01 27.63 1.00 2015 2056
## y0[7] 22.35 22.95 32.23 1.46 17.76 28.59 1.00 1712 1890
## y0[8] 18.04 18.27 9.61 1.33 12.62 23.38 1.00 2112 1851
## y0[9] 17.55 15.07 68.39 1.24 10.32 20.18 1.00 1947 1964
## y0[10] 5.00 6.17 47.07 1.40 0.78 11.56 1.00 1921 1973
## y0[11] 14.72 14.64 16.08 1.27 9.16 19.42 1.00 1831 1759
## y0[12] 16.68 16.90 14.44 1.28 11.67 21.51 1.00 1969 1988
## y0[13] 20.47 20.61 23.53 1.39 15.16 25.88 1.00 1972 1861
## y0[14] 2.42 7.34 196.97 1.33 2.66 13.42 1.00 1841 1963
## y0[15] 13.15 11.66 53.16 1.18 7.51 16.48 1.00 1700 1955
## y0[16] 15.45 15.37 12.25 1.27 10.70 20.82 1.00 1812 1940
## y0[17] 19.08 19.90 42.43 1.33 14.72 24.61 1.00 1877 1928
## y0[18] 23.19 22.19 50.14 1.43 17.12 27.88 1.00 2255 1870
## y0[19] 7.56 7.54 20.99 1.43 2.53 13.63 1.00 1768 1748
## y0[20] 18.11 21.30 140.14 1.35 15.61 27.39 1.00 2027 1974
## m1[1] 6.16 6.16 0.53 0.40 5.29 7.00 1.01 487 349
## m1[2] 8.05 8.05 0.45 0.35 7.32 8.76 1.01 499 351
## m1[3] 9.94 9.95 0.38 0.30 9.33 10.52 1.01 537 406
## m1[4] 11.82 11.84 0.33 0.28 11.28 12.31 1.01 651 456
## m1[5] 13.71 13.73 0.30 0.27 13.21 14.15 1.01 976 586
## m1[6] 15.60 15.62 0.30 0.30 15.07 16.05 1.00 1583 959
## m1[7] 17.49 17.50 0.34 0.34 16.91 18.01 1.00 2140 1465
## m1[8] 19.37 19.39 0.40 0.41 18.70 19.99 1.00 2041 1645
## m1[9] 21.26 21.28 0.47 0.48 20.47 22.00 1.00 1712 1621
## m1[10] 23.15 23.16 0.55 0.57 22.24 24.01 1.00 1407 1387
## y1[1] 6.82 6.19 55.80 1.37 0.40 11.22 1.00 1635 1644
## y1[2] 8.06 8.03 30.11 1.26 2.16 13.63 1.00 1891 1767
## y1[3] 8.79 9.93 29.38 1.31 4.76 15.37 1.00 1936 1890
## y1[4] 11.18 11.82 22.64 1.25 6.15 16.99 1.00 1856 1940
## y1[5] 11.07 13.72 181.16 1.25 9.01 18.90 1.00 1791 1713
## y1[6] 15.33 15.63 34.36 1.27 9.20 21.03 1.00 2103 1872
## y1[7] 14.07 17.53 187.16 1.28 11.70 22.31 1.00 1840 1882
## y1[8] 19.89 19.39 25.94 1.30 13.31 25.16 1.00 2141 1846
## y1[9] 20.55 21.24 32.17 1.36 15.00 26.65 1.00 1929 1763
## y1[10] 24.10 23.08 42.41 1.44 17.58 28.46 1.00 1973 1972
censored data
y i=1-N, y~N(m,s)
acutualy ya i=1-Na
lower censored yl i=1-Nl, y<L, P(y<L)=cdf(L-m /s)
upper censored yu i=1-Nu, y>U, P(y>U)=ccdf(U-m /s)
cdf(z) cumulative normal density function, P((-Infinity,z],z~N(0,1))
ccdf(z) complementary CDF, P([z,Infinity),z~N(0,1))
P(y | x,m,s)=P(ya i=1-Na)* P(yl i=1-Nl)* P(yu i=1-Nu)ex11.stan
objective variable is censored
data{
int N;
vector[N] x;
vector[N] y;
real L;
int Nl;
vector[Nl] xl;
real U;
int Nu;
vector[Nu] xu;
int N1;
vector[N1] x1;
}
parameters{
real b0;
real b1;
real<lower=0> s;
}
model{
y~normal(b0+b1*x,s);
for(i in 1:Nl)
target+=normal_lcdf(L | b0+b1*xl[i],s);
for(i in 1:Nu)
target+=normal_lccdf(U | b0+b1*xu[i],s);
}
generated quantities{
vector[N] m0;
vector[N] y0;
for(i in 1:N){
m0[i]=b0+b1*x[i];
y0[i]=normal_rng(m0[i],s);
}
vector[N1] m1;
vector[N1] y1;
for(i in 1:N1){
m1[i]=b0+b1*x1[i];
y1[i]=normal_rng(m1[i],s);
}
}n0=20
x=runif(n0,0,9)
y=rnorm(n0,10+3*x,4)
L=15
y[y<L]=L
nl=sum(y==L)
U=30
y[y>U]=U
nu=sum(y==U)
n=n0-nl-nu
qplot(x,y)
xy0=tibble(x=x,y=y)
xya=filter(xy0, y>L & y<U)
xyl=filter(xy0, y==L)
xyu=filter(xy0, y==U)
n1=10
x1=seq(min(x),max(x),length.out=n1) # new data fpr predcitiondata=list(N=n,x=xya$x,y=xya$y,N1=n1,x1=x1)
mdl=cmdstan_model('./ex8-0.stan')
mle=mdl$optimize(data=data)## Initial log joint probability = -1250.55
## Iter log prob ||dx|| ||grad|| alpha alpha0 # evals Notes
## 26 -22.0996 0.000230068 0.000160309 1 1 34
## Optimization terminated normally:
## Convergence detected: relative gradient magnitude is below tolerance
## Finished in 0.1 seconds.mle## variable estimate
## lp__ -22.10
## b0 16.93
## b1 1.47
## s 3.32
## m0[1] 23.18
## m0[2] 22.69
## m0[3] 18.66
## m0[4] 23.58
## m0[5] 23.65
## m0[6] 19.29
## m0[7] 26.20
## m0[8] 21.62
## m0[9] 28.05
## m0[10] 28.05
## m0[11] 26.85
## m0[12] 26.94
## m0[13] 27.82
## y0[1] 25.39
## y0[2] 23.03
## y0[3] 18.22
## y0[4] 21.44
## y0[5] 27.01
## y0[6] 16.77
## y0[7] 22.29
## y0[8] 20.38
## y0[9] 27.40
## y0[10] 28.59
## y0[11] 32.01
## y0[12] 32.23
## y0[13] 26.23
## m1[1] 17.64
## m1[2] 18.98
## m1[3] 20.33
## m1[4] 21.68
## m1[5] 23.03
## m1[6] 24.38
## m1[7] 25.72
## m1[8] 27.07
## m1[9] 28.42
## m1[10] 29.77
## y1[1] 17.86
## y1[2] 21.30
## y1[3] 24.32
## y1[4] 21.85
## y1[5] 19.48
## y1[6] 16.98
## y1[7] 24.80
## y1[8] 24.39
## y1[9] 31.10
## y1[10] 36.83mcmc=goMCMC(mdl,data)## Running MCMC with 4 parallel chains...
##
## Chain 1 WARNING: There aren't enough warmup iterations to fit the
## Chain 1 three stages of adaptation as currently configured.
## Chain 1 Reducing each adaptation stage to 15%/75%/10% of
## Chain 1 the given number of warmup iterations:
## Chain 1 init_buffer = 15
## Chain 1 adapt_window = 75
## Chain 1 term_buffer = 10
## Chain 1 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 1 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 1 Iteration: 600 / 600 [100%] (Sampling)
## Chain 2 WARNING: There aren't enough warmup iterations to fit the
## Chain 2 three stages of adaptation as currently configured.
## Chain 2 Reducing each adaptation stage to 15%/75%/10% of
## Chain 2 the given number of warmup iterations:
## Chain 2 init_buffer = 15
## Chain 2 adapt_window = 75
## Chain 2 term_buffer = 10
## Chain 2 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 2 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 2 Iteration: 600 / 600 [100%] (Sampling)
## Chain 3 WARNING: There aren't enough warmup iterations to fit the
## Chain 3 three stages of adaptation as currently configured.
## Chain 3 Reducing each adaptation stage to 15%/75%/10% of
## Chain 3 the given number of warmup iterations:
## Chain 3 init_buffer = 15
## Chain 3 adapt_window = 75
## Chain 3 term_buffer = 10
## Chain 3 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 3 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 3 Iteration: 600 / 600 [100%] (Sampling)
## Chain 4 WARNING: There aren't enough warmup iterations to fit the
## Chain 4 three stages of adaptation as currently configured.
## Chain 4 Reducing each adaptation stage to 15%/75%/10% of
## Chain 4 the given number of warmup iterations:
## Chain 4 init_buffer = 15
## Chain 4 adapt_window = 75
## Chain 4 term_buffer = 10
## Chain 4 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 4 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 4 Iteration: 600 / 600 [100%] (Sampling)
## Chain 1 finished in 0.0 seconds.
## Chain 2 finished in 0.0 seconds.
## Chain 3 finished in 0.0 seconds.
## Chain 4 finished in 0.0 seconds.
##
## All 4 chains finished successfully.
## Mean chain execution time: 0.0 seconds.
## Total execution time: 0.2 seconds.seeMCMC(mcmc,'m')## variable mean median sd mad q5 q95 rhat ess_bulk ess_tail
## lp__ -22.70 -22.27 1.65 1.15 -25.85 -21.09 1.02 300 354
## b0 16.74 16.72 3.33 2.88 11.88 21.37 1.03 263 333
## b1 1.50 1.51 0.61 0.53 0.60 2.38 1.02 289 390
## s 4.14 3.94 1.13 0.87 2.80 6.16 1.01 757 491
## m0[1] 23.11 23.10 1.32 1.19 20.96 25.17 1.01 580 634
## m0[2] 22.61 22.61 1.40 1.25 20.30 24.76 1.02 464 530
## m0[3] 18.50 18.47 2.67 2.27 14.50 22.21 1.03 269 331
## m0[4] 23.51 23.52 1.26 1.14 21.48 25.49 1.01 722 893
## m0[5] 23.58 23.58 1.26 1.13 21.57 25.55 1.01 753 892
## m0[6] 19.14 19.13 2.44 2.07 15.44 22.57 1.03 273 334
## m0[7] 26.18 26.23 1.45 1.29 23.81 28.47 1.00 1399 844
## m0[8] 21.52 21.53 1.67 1.46 18.88 24.02 1.02 337 326
## m0[9] 28.06 28.10 1.97 1.74 24.81 31.15 1.01 598 554
## m0[10] 28.06 28.10 1.97 1.74 24.81 31.15 1.01 599 554
## m0[11] 26.84 26.90 1.61 1.39 24.23 29.34 1.00 901 704
## m0[12] 26.94 27.00 1.63 1.43 24.27 29.49 1.00 857 700
## m0[13] 27.83 27.87 1.90 1.68 24.67 30.77 1.01 633 586
## y0[1] 23.25 23.24 4.49 4.13 15.98 30.68 1.00 1854 1717
## y0[2] 22.61 22.69 4.59 4.17 15.13 29.71 1.00 1688 1365
## y0[3] 18.48 18.60 4.99 4.36 10.57 26.37 1.01 798 885
## y0[4] 23.47 23.50 4.56 4.22 16.12 30.55 1.00 1836 1799
## y0[5] 23.62 23.62 4.44 4.16 16.49 30.55 1.00 1767 1702
## y0[6] 19.26 19.34 4.79 4.48 11.47 26.62 1.01 898 1250
## y0[7] 26.34 26.32 4.52 4.14 19.10 33.59 1.00 1778 1533
## y0[8] 21.51 21.62 4.66 4.27 13.41 28.95 1.00 1254 1364
## y0[9] 27.99 27.96 4.93 4.35 20.15 35.73 1.00 1604 1262
## y0[10] 28.06 28.07 4.84 4.41 20.36 35.22 1.00 1451 1271
## y0[11] 26.64 26.67 4.59 4.15 19.08 34.33 1.00 1850 1891
## y0[12] 26.98 27.14 4.88 4.24 19.14 34.82 1.00 1842 1593
## y0[13] 27.82 27.86 4.67 4.16 20.50 35.36 1.00 1415 1186
## m1[1] 17.46 17.43 3.05 2.64 12.98 21.70 1.03 264 316
## m1[2] 18.83 18.81 2.55 2.16 14.97 22.38 1.03 270 326
## m1[3] 20.20 20.21 2.07 1.80 17.00 23.18 1.02 290 307
## m1[4] 21.58 21.60 1.65 1.44 18.97 24.06 1.02 341 329
## m1[5] 22.95 22.95 1.34 1.20 20.77 25.02 1.01 541 614
## m1[6] 24.32 24.34 1.22 1.10 22.35 26.21 1.00 1216 978
## m1[7] 25.69 25.74 1.35 1.24 23.50 27.86 1.00 1674 928
## m1[8] 27.07 27.12 1.67 1.45 24.32 29.67 1.00 809 710
## m1[9] 28.44 28.48 2.09 1.83 25.01 31.68 1.01 554 575
## m1[10] 29.81 29.83 2.57 2.28 25.71 33.77 1.01 457 582
## y1[1] 17.46 17.46 5.37 4.65 9.06 26.00 1.01 629 771
## y1[2] 18.74 18.77 4.99 4.50 10.96 26.75 1.01 823 804
## y1[3] 20.08 20.18 4.81 4.51 12.50 27.73 1.00 1135 1265
## y1[4] 21.58 21.63 4.77 4.14 13.90 29.31 1.01 1047 1630
## y1[5] 23.08 23.05 4.55 4.37 16.10 30.64 1.00 1449 1428
## y1[6] 24.30 24.21 4.31 4.06 17.62 31.16 1.00 1820 1738
## y1[7] 25.69 25.58 4.50 3.98 18.48 33.11 1.00 2025 1510
## y1[8] 27.08 27.17 4.53 4.16 19.58 34.17 1.00 1703 1606
## y1[9] 28.30 28.27 4.70 4.27 20.56 35.96 1.00 1780 1776
## y1[10] 29.92 29.94 5.15 4.59 21.62 38.20 1.01 1044 1167
y0=mcmc$draws('y0')
smy0=summary(y0)
grid.arrange(
qplot(xya$y,smy0$median,xlab='obs.',ylab='prd.')+geom_abline(intercept=0,slope=1))
par(mfrow=c(1,2))
hist(xya$y-smy0$median,xlab='obs.-prd.',main='residual')
density(xya$y-smy0$median) |> plot(xlab='obs.-prd.',main='residual')
m1=mcmc$draws('m1')
smm1=summary(m1)
y1=mcmc$draws('y1')
smy1=summary(y1)
xy=tibble(x=x1,m=smm1$median,ml5=smm1$q5,mu5=smm1$q95,yl5=smy1$q5,yu5=smy1$q95)
qplot(x,y,col=I('red'))+
geom_line(aes(x=x,y=ml5),xy,col='darkgray')+
geom_line(aes(x=x,y=mu5),xy,col='darkgray')+
geom_line(aes(x=x,y=yl5),xy,col='lightgray')+
geom_line(aes(x=x,y=yu5),xy,col='lightgray')+
geom_line(aes(x=x,y=m),xy,col='black')
data=list(N=n,x=xya$x,y=xya$y,
L=L,Nl=nl,xl=xyl$x,
U=U,Nu=nu,xu=xyu$x,
N1=n1,x1=x1)
mdl=cmdstan_model('./ex11.stan')
mle=mdl$optimize(data=data)## Rejecting initial value:
## Log probability evaluates to log(0), i.e. negative infinity.
## Stan can't start sampling from this initial value.
## Rejecting initial value:
## Log probability evaluates to log(0), i.e. negative infinity.
## Stan can't start sampling from this initial value.
## Rejecting initial value:
## Log probability evaluates to log(0), i.e. negative infinity.
## Stan can't start sampling from this initial value.
## Rejecting initial value:
## Log probability evaluates to log(0), i.e. negative infinity.
## Stan can't start sampling from this initial value.
## Rejecting initial value:
## Log probability evaluates to log(0), i.e. negative infinity.
## Stan can't start sampling from this initial value.
## Rejecting initial value:
## Log probability evaluates to log(0), i.e. negative infinity.
## Stan can't start sampling from this initial value.
## Rejecting initial value:
## Log probability evaluates to log(0), i.e. negative infinity.
## Stan can't start sampling from this initial value.
## Rejecting initial value:
## Log probability evaluates to log(0), i.e. negative infinity.
## Stan can't start sampling from this initial value.
## Rejecting initial value:
## Log probability evaluates to log(0), i.e. negative infinity.
## Stan can't start sampling from this initial value.
## Initial log joint probability = -99.1728
## Iter log prob ||dx|| ||grad|| alpha alpha0 # evals Notes
## 26 -28.7063 0.000276867 2.73144e-05 1 1 28
## Optimization terminated normally:
## Convergence detected: relative gradient magnitude is below tolerance
## Finished in 0.1 seconds.mle## variable estimate
## lp__ -28.71
## b0 11.13
## b1 2.46
## s 4.08
## m0[1] 21.57
## m0[2] 20.76
## m0[3] 14.02
## m0[4] 22.23
## m0[5] 22.35
## m0[6] 15.07
## m0[7] 26.61
## m0[8] 18.96
## m0[9] 29.69
## m0[10] 29.69
## m0[11] 27.69
## m0[12] 27.85
## m0[13] 29.32
## y0[1] 26.99
## y0[2] 22.04
## y0[3] 16.42
## y0[4] 21.22
## y0[5] 20.17
## y0[6] 17.13
## y0[7] 27.84
## y0[8] 20.35
## y0[9] 23.36
## y0[10] 32.12
## y0[11] 26.68
## y0[12] 25.98
## y0[13] 31.74
## m1[1] 12.31
## m1[2] 14.56
## m1[3] 16.81
## m1[4] 19.06
## m1[5] 21.31
## m1[6] 23.56
## m1[7] 25.81
## m1[8] 28.06
## m1[9] 30.31
## m1[10] 32.56
## y1[1] 11.28
## y1[2] 10.30
## y1[3] 17.24
## y1[4] 20.66
## y1[5] 18.76
## y1[6] 14.29
## y1[7] 25.87
## y1[8] 26.28
## y1[9] 23.49
## y1[10] 32.09mcmc=goMCMC(mdl,data)## Running MCMC with 4 parallel chains...
##
## Chain 1 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 1 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 1 Iteration: 600 / 600 [100%] (Sampling)
## Chain 2 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 2 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 2 Iteration: 600 / 600 [100%] (Sampling)
## Chain 3 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 3 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 3 Iteration: 600 / 600 [100%] (Sampling)
## Chain 4 WARNING: There aren't enough warmup iterations to fit the
## Chain 4 three stages of adaptation as currently configured.
## Chain 4 Reducing each adaptation stage to 15%/75%/10% of
## Chain 4 the given number of warmup iterations:
## Chain 4 init_buffer = 15
## Chain 4 adapt_window = 75
## Chain 4 term_buffer = 10
## Chain 4 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 4 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 4 Iteration: 600 / 600 [100%] (Sampling)
## Chain 1 finished in 0.1 seconds.
## Chain 2 finished in 0.1 seconds.
## Chain 3 finished in 0.0 seconds.
## Chain 4 finished in 0.0 seconds.
##
## All 4 chains finished successfully.
## Mean chain execution time: 0.1 seconds.
## Total execution time: 0.2 seconds.seeMCMC(mcmc,'m',ch=T)## variable mean median sd mad q5 q95 rhat ess_bulk ess_tail
## lp__ -29.03 -28.64 1.44 1.14 -32.08 -27.49 1.01 424 553
## b0 10.10 10.41 3.10 2.76 4.50 14.27 1.02 237 253
## b1 2.65 2.58 0.57 0.51 1.87 3.73 1.01 293 258
## s 5.16 4.87 1.41 1.22 3.43 7.76 1.00 813 932
## m0[1] 21.36 21.45 1.37 1.22 18.97 23.43 1.00 504 588
## m0[2] 20.48 20.59 1.44 1.28 17.89 22.59 1.00 405 457
## m0[3] 13.22 13.50 2.51 2.22 8.70 16.62 1.02 239 260
## m0[4] 22.07 22.15 1.33 1.22 19.77 24.09 1.00 629 692
## m0[5] 22.20 22.26 1.33 1.21 19.90 24.20 1.00 655 719
## m0[6] 14.35 14.59 2.31 2.07 10.16 17.53 1.01 242 274
## m0[7] 26.80 26.78 1.51 1.38 24.39 29.37 1.01 1473 1131
## m0[8] 18.55 18.72 1.65 1.45 15.49 20.92 1.01 300 331
## m0[9] 30.12 30.03 1.96 1.83 27.11 33.54 1.01 905 909
## m0[10] 30.12 30.02 1.96 1.83 27.11 33.54 1.01 906 909
## m0[11] 27.96 27.92 1.65 1.53 25.34 30.79 1.01 1275 990
## m0[12] 28.14 28.09 1.67 1.55 25.46 31.03 1.01 1240 966
## m0[13] 29.72 29.63 1.90 1.79 26.80 33.05 1.01 968 875
## y0[1] 21.42 21.59 5.43 5.00 12.48 30.22 1.00 2006 1586
## y0[2] 20.28 20.46 5.62 5.29 10.93 29.15 1.00 1928 1778
## y0[3] 13.20 13.34 5.96 5.51 3.32 22.70 1.00 1250 1329
## y0[4] 21.94 22.05 5.35 4.91 13.43 30.42 1.00 1822 1973
## y0[5] 22.48 22.56 5.77 5.08 13.32 31.31 1.00 2020 1431
## y0[6] 14.52 14.80 5.89 5.35 4.23 23.33 1.00 1109 1231
## y0[7] 26.77 26.72 5.46 5.02 17.84 35.73 1.00 1983 1913
## y0[8] 18.43 18.53 5.60 5.01 9.24 27.24 1.00 1394 1568
## y0[9] 30.10 30.03 5.62 5.01 21.27 39.47 1.00 1813 1895
## y0[10] 30.07 29.96 5.84 5.19 20.34 39.64 1.00 1797 1754
## y0[11] 27.81 27.72 5.77 5.09 18.77 37.37 1.00 1792 1750
## y0[12] 28.24 28.27 5.57 5.18 19.23 37.26 1.00 2034 1698
## y0[13] 29.70 29.54 5.67 5.26 20.70 39.36 1.00 1886 1242
## m1[1] 11.37 11.67 2.86 2.53 6.23 15.22 1.02 237 261
## m1[2] 13.80 14.06 2.41 2.14 9.46 17.08 1.01 240 273
## m1[3] 16.23 16.45 1.99 1.79 12.52 19.01 1.01 255 269
## m1[4] 18.65 18.83 1.64 1.43 15.62 21.02 1.01 303 342
## m1[5] 21.08 21.17 1.39 1.23 18.64 23.14 1.00 465 588
## m1[6] 23.51 23.51 1.31 1.21 21.35 25.46 1.01 1032 990
## m1[7] 25.94 25.91 1.43 1.29 23.67 28.34 1.01 1509 1139
## m1[8] 28.36 28.31 1.70 1.59 25.67 31.33 1.01 1197 976
## m1[9] 30.79 30.67 2.07 1.91 27.61 34.46 1.01 819 904
## m1[10] 33.22 33.02 2.50 2.34 29.51 37.76 1.01 629 601
## y1[1] 11.38 11.77 6.02 5.29 0.89 20.60 1.00 864 1156
## y1[2] 13.97 14.23 6.02 5.51 3.79 23.28 1.00 1089 1142
## y1[3] 16.17 16.33 5.81 5.12 6.35 25.35 1.00 1174 1172
## y1[4] 18.58 18.70 5.73 5.18 9.36 27.49 1.00 1207 1298
## y1[5] 21.20 21.20 5.57 4.99 11.90 30.27 1.00 1685 1654
## y1[6] 23.52 23.55 5.51 5.01 14.48 32.31 1.00 1902 1961
## y1[7] 25.94 25.82 5.61 5.26 16.85 35.08 1.00 1895 1638
## y1[8] 28.34 28.30 5.56 5.08 19.67 37.42 1.00 1760 1434
## y1[9] 30.83 30.67 5.62 5.16 21.97 39.99 1.00 1297 1379
## y1[10] 33.20 33.23 5.86 5.24 23.78 42.71 1.00 1725 1782
y0=mcmc$draws('y0')
smy0=summary(y0)
grid.arrange(
qplot(xya$y,smy0$median,xlab='obs.',ylab='prd.')+geom_abline(intercept=0,slope=1))
par(mfrow=c(1,2))
hist(xya$y-smy0$median,xlab='obs.-prd.',main='residual')
density(xya$y-smy0$median) |> plot(xlab='obs.-prd.',main='residual')
m1=mcmc$draws('m1')
smm1=summary(m1)
y1=mcmc$draws('y1')
smy1=summary(y1)
xy=tibble(x=x1,m=smm1$median,ml5=smm1$q5,mu5=smm1$q95,yl5=smy1$q5,yu5=smy1$q95)
qplot(x,y,col=I('red'))+
geom_line(aes(x=x,y=ml5),xy,col='darkgray')+
geom_line(aes(x=x,y=mu5),xy,col='darkgray')+
geom_line(aes(x=x,y=yl5),xy,col='lightgray')+
geom_line(aes(x=x,y=yu5),xy,col='lightgray')+
geom_line(aes(x=x,y=m),xy,col='black')
sensitivity/specificity analysis
ex14.stan
estimating sens and spec
data {
int N;
array[N] int x;
array[N] int y;
}
parameters {
real<lower=0,upper=1> p;
real<lower=0,upper=1> se;
real<lower=0,upper=1> sp;
}
model {
p~uniform(0,1);
se~uniform(0,1);
sp~uniform(0,1);
for (i in 1:N) {
y[i]~bernoulli(x[i]*se+(1-x[i])*(1-sp));
}
}
generated quantities {
real ppv;
real npv;
ppv=se*p/((se*p)+((1-p)*(1-sp)));
npv=(1-p)*sp/(((1-p)*sp)+(p*(1-se)));
}n=20
data=list(N=n,
x=sample(0:1,n,replace=T),
y=sample(0:1,n,replace=T))
mdl=cmdstan_model('./ex14.stan')
mle=mdl$optimize(data=data)## Initial log joint probability = -16.7538
## Iter log prob ||dx|| ||grad|| alpha alpha0 # evals Notes
## 5 -13.7528 0.00030773 2.0269e-05 1 1 8
## Optimization terminated normally:
## Convergence detected: relative gradient magnitude is below tolerance
## Finished in 0.1 seconds.mle## variable estimate
## lp__ -13.75
## p 0.78
## se 0.54
## sp 0.43
## ppv 0.77
## npv 0.21system.time({
mcmc=goMCMC(mdl,data)
seeMCMC(mcmc,ch=F)
})## Running MCMC with 4 parallel chains...
##
## Chain 1 WARNING: There aren't enough warmup iterations to fit the
## Chain 1 three stages of adaptation as currently configured.
## Chain 1 Reducing each adaptation stage to 15%/75%/10% of
## Chain 1 the given number of warmup iterations:
## Chain 1 init_buffer = 15
## Chain 1 adapt_window = 75
## Chain 1 term_buffer = 10
## Chain 1 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 1 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 1 Iteration: 600 / 600 [100%] (Sampling)
## Chain 2 WARNING: There aren't enough warmup iterations to fit the
## Chain 2 three stages of adaptation as currently configured.
## Chain 2 Reducing each adaptation stage to 15%/75%/10% of
## Chain 2 the given number of warmup iterations:
## Chain 2 init_buffer = 15
## Chain 2 adapt_window = 75
## Chain 2 term_buffer = 10
## Chain 2 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 2 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 2 Iteration: 600 / 600 [100%] (Sampling)
## Chain 3 WARNING: There aren't enough warmup iterations to fit the
## Chain 3 three stages of adaptation as currently configured.
## Chain 3 Reducing each adaptation stage to 15%/75%/10% of
## Chain 3 the given number of warmup iterations:
## Chain 3 init_buffer = 15
## Chain 3 adapt_window = 75
## Chain 3 term_buffer = 10
## Chain 3 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 3 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 3 Iteration: 600 / 600 [100%] (Sampling)
## Chain 4 WARNING: There aren't enough warmup iterations to fit the
## Chain 4 three stages of adaptation as currently configured.
## Chain 4 Reducing each adaptation stage to 15%/75%/10% of
## Chain 4 the given number of warmup iterations:
## Chain 4 init_buffer = 15
## Chain 4 adapt_window = 75
## Chain 4 term_buffer = 10
## Chain 4 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 4 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 4 Iteration: 600 / 600 [100%] (Sampling)
## Chain 1 finished in 0.0 seconds.
## Chain 2 finished in 0.0 seconds.
## Chain 3 finished in 0.0 seconds.
## Chain 4 finished in 0.0 seconds.
##
## All 4 chains finished successfully.
## Mean chain execution time: 0.0 seconds.
## Total execution time: 0.2 seconds.
##
## variable mean median sd mad q5 q95 rhat ess_bulk ess_tail
## lp__ -19.54 -19.20 1.29 1.01 -22.07 -18.13 1.00 718 912
## p 0.51 0.50 0.29 0.37 0.05 0.95 1.01 609 458
## se 0.53 0.53 0.12 0.12 0.34 0.72 1.00 3217 1505
## sp 0.45 0.44 0.16 0.17 0.20 0.71 1.00 2121 1464
## ppv 0.50 0.50 0.30 0.39 0.05 0.95 1.01 632 479
## npv 0.48 0.47 0.29 0.38 0.05 0.95 1.01 633 539
## ユーザ システム 経過
## 0.558 0.121 0.699appendix: bimodal distribution, mixed normal distribution
stan
data {
int<lower=0> N;
real y[N];
}
parameters {
real<lower=0, upper=1> theta; //ratio of mix
real mu1;
real mu2;
real<lower=0> sigma1;
real<lower=0> sigma2;
}
model {
for (n in 1:N) {
target += log_mix(theta,
normal_lpdf(y[n] | mu1, sigma1),
normal_lpdf(y[n] | mu2, sigma2));
}
}EM algorithm
y0=rnorm(100,0,1)
y1=rnorm(100,-5,1)
y2=rnorm(100,5,1)
y=sample(c(y0,y1,y2),100)
density(y) |> plot()
library(mclust)
rst=Mclust(y)
summary(rst)## ----------------------------------------------------
## Gaussian finite mixture model fitted by EM algorithm
## ----------------------------------------------------
##
## Mclust E (univariate, equal variance) model with 3 components:
##
## log-likelihood n df BIC ICL
## -239 100 6 -506 -507
##
## Clustering table:
## 1 2 3
## 28 44 28rst$parameters## $pro
## [1] 0.281 0.438 0.281
##
## $mean
## 1 2 3
## -5.214 -0.047 5.029
##
## $variance
## $variance$modelName
## [1] "E"
##
## $variance$d
## [1] 1
##
## $variance$G
## [1] 3
##
## $variance$sigmasq
## [1] 0.842
##
##
## $Vinv
## NULLplot(rst)
